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I. INTRODUCTION 


Thermoacoustic heat transport is a process through which 
an acoustic field generates, or inversely is generated fron, 
a flow of heat. Like every engine, a thermoacoustic engine can 
be configured as either type of classical heat engine, a heat 
pump (or refrigerator) or a prime mover. The purpose of this 
thesis is the computer implementation of a general formulation 
of thermoacoustics developed by Arnott et al [Ref. 1] using 
Matlab. In particular, we find the theoretical values of the 
quality factor and the resonance frequency of a prime mover as 
a function of the temperature difference applied across the 
prime mover stack. The program determines the impedance and 
normalized pressure distribution along the prime mover 
beginning from a known value of specific acoustic impedance 
and normalized pressure amplitude at one end. It finds the 
resonance frequency and quality factor through calculation of 
the pressure amplitude at the opposite end as a function of 
frequency. The result of the computer implementation are 
compared with measurements made with open and closed ended 
prime mover [Ref. 2,3], and with the results of a standing 


wave analysis of prime movers [Ref. 4]. 


II. THEORY 
The goal of this chapter is to describe the basic 
geometry of a prime mover and to summarize the important 
points of Arnott’s method. Parallel with this we explain how 
this theory is implemented in the computer program. The reader 


is referred to [{Ref. 1] for full details of Arnott’s method. 


A. ARNOTT’S METHOD. 

The necessary background for this thesis can be explained 
with the help of Fig.l which shows the basic constructioneee 
the thermoacoustic prime mover. It constists of a circular 
cross section acoustic resonator , rigidly capped at both 
ends. Situated within the resonator are the thermoacoustic 
elements consisting of an ambient heat exchanger, a prime 
mover stack anda hot heat exchanger. In our prime mover both 
the stack and heat exchangers are parallel plates spaced by a 
few thermal penetration depths. Complete details of the prime 
mover construction are given in [Ref. 2,3,4]. The traverse 
coordinates of the prime mover are taken to be x and y, and 
the longitudinal coordinate is z. The hot and cold sections 
are Circular ducts open at the heat exchanger end and open or 
closed at the other end. The prime mover is divided into 
either 5 or 7 sections depending on the boundary conditions 


(Fig. 1,2]. The stack is divided into 11 subsections. 
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Figure 1. Rigid end Prime Mover 
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Figure 2. Open end Prime Mover 


The computer program finds the normalized acoustic 
pressure and the specific acoustic impedance in all parts of 
the prime mover as a function of frequency. The pressure as a 
function of frequency gives the frequency response. From the 
frequency response we can determine the quality factor (Q) and 
the resonance frequency (f,). To do that we assume a linear 
temperature gradient across the stack via the hot and ambient 
(cold) heat exchangers. The stack is divided into eleven 
Subsections. A constant temperature 1S assumed in each 
subsection of the stack. In all others parts of the prime 
mover we assume that the temperature iS constant and equal to 
that of the nearest heat exchanger. The ambient temperature 
(T,) is taken to be a function of z (T,(z)) only in the stack 
area. 

The program starts with specific acoustic impedance at 
the right end. There are two boundary conditions of interest: 
rigid and open. The specific acoustic impedance of a rigid end 


ns 


=(14i) ¥2__1 Bkoce (1) 
BEA 42) Ty ay? ON wn VNpr: 


momeescribed in [Ref. 5:p.529]. In the above equation N,, is 
Pie Prandtl] number, c, is the isobaric heat capacity per unit 


mass, Y is the ratio c,/c, where c, is the constant volume heat 


capacity per unit mass, N 1s the viscosity and p, represents 
the ambient value of the density. The specific acoustic 


impedance of an open end is: 


Z=cp ka(~2-0.64), (2) 


where a is the radius of the tube and k=@/c is the wavenumber, 
® is the angular frequency and i=(-1)* [Ref. 6:p.202]. 

The acoustic pressure iS specified at the same end. 
Because the Q is independent of the absolute value of the 
acoustic pressure in linear theory, the specified pressure is 
only a relative or normalized pressure. Next we find the 
pressure and the impedance at the entrance to the hot heat 
exchanger using the Rayleigh’s impedance translation 
theorem [(Ref.1,5]. This is one difference between Arnott’s and 
other methods. That is, other methods find the fluid variables 
by integrating the governing equations along the prime mover. 
The impedance translation theorem state that: within any 
homogeneous layer, with intrinsic (or characteristic) specific 
acoustic impedance Z,,,, the local specific impedance Z(z-1l) at 


zZ-l 1s related to that at 2 by [Ret oes 


Z(y) COSIWURT) [1 2777s Taek) 


2AZ-1) =2ine7 Gos (kl) -dZly) Sin(kl) 


Likewise, the acoustic pressure at the entrance to the hot 
heat exchanger is found through the pressure translation 


theorem: 


Z; 


P, (z-1) =P, (z) {cos (k1) -i[ mi 





ne j}sin(kl)}, (4) 
Z) 


where Z(z) is the specific impedance at z, P,(z) is the first 
Peeeers Value of pressure, 1.¢€, the acoustic pressure, and 2,,, 


2S: 


Exon (5) 


Zine” TFN DK] 


In the above equation, 2 represents the porosity and is equal 
to NA, where A is the cross sectional area of single pore and 
N is the number of pores per unit area. k is the complex wave 


number in the pore and given by: 





2a-w* 1 vee 6 
(KN Wg) #222 A (y ~ Cy -2) F(A 9). (6) 


The F(A) and F(A,) factors arise from the equation for 


the z component of the acoustic velocity: 


F(x,y;A ) GP, (Zz) (7) 


U (x,y,z) = iop dz 


The F(x,y;A) satisfies the following differential equation 


[Refi asia 


2 


F(x, yj ) +( Ve (Cx) y 7a) ale (8) 





id ? 
where the Laplacian operator V’, is given by: 


v= §- §? 


+ 


8 x? 8 y2_ a 








In our case where the pores of the stack are parallel plates, 
the F depends only on y and the shear wave number A or the 


dimensionless thermal disturbance number A, which are given 





Dye 
1 
N= Rien ya (10) 
n 
and 
a 
h aA NZ. a! 
The F(y;A) is given by: 
cosh (y=i4. ¥) 
Py Xk = oo (gy 


cosh (¥=1*_) 


In the above equation a is the half separation distance 
between the parallel plates of the stack or heat exchanger. 
The pore average of F(y;A) for heat exchangers and stack, is 


given by: 


F(A )=1-(2 y=a) tanh (/=7%.). (13) 


Finally we use the boundary layer approximation in the 


Circular regions of the prime mover: 


lim, F(A )=1-(8 /R) (1+4), (14) 


In the above relation R is twice the hydraulic radius of the 
pore, or twice the ratio of the transverse pore area to the 
pore perimeter. Actually for heat exchangers and stack with 
parallel plates, Ris the separation distance (2a) between the 
plates and for the circular parts it is the radius of the 
tube. 8 is the viscous or thermal penetration depth. 

emcee FP, and Z are known at the entrance to heat exchanger, 
we need the values just inside the heat exchanger. Using the 
pressure and the impedance translation theorem we have a 
solution up to the stack where a linear temperature gradients 
exists. In the heat exchanger three things are different from 


the circular parts. These are: the wavenumbers in the heat 


exchangers are complex, the porosity now is NA, and F(A) as we 
mention above is given by equation (13). 

The point now is how to find pressure and impedance in 
the stack. In this case Arnott’s method uses the following 
system of the first order differential equations that can be 


solved using Range-Kutta methods: 


azZitale: Ez) * 
az 1K (2) Zine (2) ca (Zea (15) 
and 
Gia ae CZs) 16 
TaiGioInn ne eo eae ( ) 
where Q(z) is: 
121 a 
oe (17) 
a (A, Aq) a8 [EA] 
T.. 1s the temperature gradient given by: 
aT, (z) 
=e (18) 
To2 dz 


8B is the thermal expansion coefficient: 
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ae 9 TP (19) 
p 


In our program, this system of equations is solved with the 
Runge-Kutta-Fehlberg method for every subsection of the stack. 
In so doing, the impedance and the pressure distribution in 
the stack is obtained. Now using the pressure and the 
impedance translation theorem we find the pressure and the 
impedance in the rest of the prime mover. 

Knowing the pressure at the left end we can calculate the 
quality factor, the resonance frequency and the maximum 
pressure amplitude. This is done by minimizing the following 


equation (Ref. 8]: 


(IPG) |fo)2 
or -amplit?, (20) 
cee Ce 
ro E, (ex 
where Q is the quality factor, |P(1)|] is the maximum 


amplitude of pressure at the left end, f, is the resonance 
frequency, amplit is the calculated amplitude of pressure at 
the left end and is a function of the frequency f. In the 
above relation there are two known variables f and amplit and 


three unknowns (Q,f,,P(1)) that we want to find. To find these 
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unknowns, the 3-parameter least Squares technique is used with 
the help of simplex method. 

Up to this point, Q and f, have been determined only by 
considering the pressure amplitude at the left end. We have 
not worried about matching the acoustic impedance at the left 
end. Impedance matching allows us to refine the values of Q 
and f£, through the complex eigenfrequency. Using this 
technique, one assumes the angular frequency is complex and 


given by: 


w =2n fy(1-) | (21) 
Substituting the initial value of f, and Q determined from the 
frequency response, into this equation gives an initial guess 
at the complex eigenfrequency for the system. This frequency 
is used to calculate the impedance throughout the prime mover 
as explained earlier. The value of impedance in the gas at the 
left end is compared to the value for the boundary condition 
(i.e. the impedance of a rigid end with thermal loss). The 
difference between the computed and actual impedance is used 
to adjust the eigenfrequency. The impedance is computed with 
the new eigenfrequency until both real and imaginary parts of 


the impedance matche at the left end. The final values of f, 


rz 


and Q are determined from the real and imaginary parts of the 


eigenfrequency that best matche the impedance. 


B. IMPLEMENTATION OF ARNOTT’S METHOD 

The program uses functions that already exist in Matlab 
where possible. For the Runge-Kutta-Fehlberg 5th order method 
we modified the existing function in Matlab to integrate with 
initial values of the integration variable greater than the 
final. This 1S necessary because we put the origin at the left 
end (z=0) and we begin our calculations from the right end 
(z=L). One reason that we use this particular Runge-Kutta 
method 1s itS accuracy. The Runge-Kutta-Fehlberg computes two 
Runge-Kutta estimates for the value y,,, but of different 
orders of error. The global error in this numerical method is 
ey and mine Local error 15 O(h°) (Ref. 7}. 

To do the minimization, we use the function fmins with 
the help of simplex numerical method. The fmin function with 
Simplex method is used to match the impedances at the left end 


using the technique of least Squares. 


ils: 


III. THE MATLAB PROGRAM 
The goal of this chapter is to provide an overview of the 
Structure of the Matlab program that implements Arnott’s 
method. Furthermore it describes all the functions, special 
tricks and limitations of this code. It also explains how to 


run the program. 


A. OUTLINE OF THE MATLAB PROGRAM 

The program is divided into three parts. In the first 
part the variables are specified for every section of the 
prime mover. They are: the type of section (open, heat 
exchanger, stack), the number of subsections that exists 
inside the section, the temperature at the left and right ends 
of the section, twice the hydraulic radius of the section, and 
the porosity of the section. Also, specified are the frequency 
range and the number of intervals in this range. 

The second part is the main program. It finds the 
impedance and the pressure distribution in every section of 
the prime mover for each frequency within the assumed range. 
To accomplish this it begins at the right end of the tube 
where it calculates the impedance of the end, rigid or open, 
and assigns an arbitrary value to the pressure. Continuing, it 
calculates using the translation theorems, the pressure and 


the impedance in all other sections except in the stack. In 
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the stack, where there exists a temperature gradient, the 
program calls the function "“trode45" (a modified version of 
the function of Matlab "“ode45") to integrate from the right 
end to the left end of the stack. This function solves 
differential equations uSing the 5th order Runge-Kutta 
numerical method with Fehlberg coefficients. 

The last part of the Matlab program uses the calculated 
pressure at the left end from the previous part to find the 
maximum amplitude, the resonance frequency and the quality 
factor. These parameters are found from a least Square fit to 
equation (20). The minimization is performed with the matlab 
Memeelon “fmins" with “options 1,2,3,14". Finally, the Q and 
f, estimates are used as a Starting points for the final 
calculation of the quality factor and the resonance frequency 
using the complex eigenfrequency method. This refinement 
usually is a small correction (in the fourth decimal place) to 
initial guesses. Finally, the program plots 1/Q and f, versus 


temperature difference and saves the results ina data file. 


B. LIMITATIONS OF THE PROGRAM, SPECIAL POINTS. 

There are two programs, one for the open end prime mover 
and one for the rigid end tube. To run the open end progam, 
peemeneeas to run “arffl.m" with the functions “error,m- 
erro3.m-argo3f.m". To run the rigid end program one needs to 


Bm argoZ2.m” with the functions “error.m-erro3.m—-argo3.m". 


iD 


The difference in the two progams is that one uses "“argo3.m" 
and the other uses "“argof3.m". Before running one program or 
the other, one needs to switch these functions in the 
"erro3.m"” E£umerlon- 

Futhermore, there are some factors that influence the 
operation of the program. One important step is the choice of 
the initial range of the frequencies. To begin the program one 
needs to know the approximate resonance frequency as well as 
which longitudinal acoustic mode is to be used. Initially the 
frequency range is f, + 20Hz. After the program runs for the 
first value of AT the frequency range is automatically 
adjusted to be: 

A f= — 
To] 


Another important factor is of course the number of 


(22) 


points between the maximum and minimum frequency that the 
program uses. Usually this number is 500 but sometimes when 
the values of quality factor gets large this may have to be 
doubled (1000 points). In these cases, it may also be 
necessary to use steps of 5 degrees Kelvin instead of the 
usual 10 degrees Kelvin. 

Finally, another important point concerns the “counters" 
that are used by the program. The most important counter 1S 
the "flag". If somebody wants to run the program for a case 
where the quality factor increases with AT, one needs to put 


the initial value of the "“flag" counter equal to zero. 
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Otherwise, it should be set equal to one. The purpose of this 
counter is to force the program to scan for negative values of 


the inverse quality factor at the proper time. 


ee? 


IV. VALIDATION AND RESULTS 
In this chapter we compare the results of the program to 
experimental results for two different prime movers [Ref. 2 
and 3} as well as with the results of the standing wave 
analysis by Atchley [Ref 4]. 

The first results are for a rigid end prime mover 
described in Ref. 2 and 4. Comparison of both 1/Q and 
resonance frequency are made. Figures 3 through 5 show 1/Q vs 
AT (in degrees K) for three different mean gas pressures. The 
gas is helium. In these figures, the open circles represent 
measured values, the *’s represent the results of the standing 
wave analysis and the solid line represents the result of the 
Matlab program. The overall agreement between the present 
analysis and the measurements is good, but not as good as for 
the standing wave analysis. The present analysis agrees best 
at low values of AT. The reasons for the increased discrepancy 
at higher values of AT are not understood. One would expect 
the present analysis to agree better with measurements than 
the standing wave analysis. This is expected because the 
present analysis makes fewer assumptions. In the particular, 
it does a better job of taking into account changes in the 
acoustic pressure amplitude in differents parts of the prime 


mover. The standing wave analysis ignores any changes in 
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pressure amplitude due to changes in cross sectional area of 
the prime mover. 

Figures 6 through 8 show results of resonance frequency 
(in Hz) vs AT for the same three mean pressures. The agreement 
between the present analysis and the measured frequencies is 
within 1% in all cases. The standing wave analysis also agrees 
with the measured values to about the same accuracy. However, 
the two methods show different depedences on AT. This is 
because 1) the standing wave analysis considers only the 
temperature dependence of the sound speed, whereas the present 
analysis includes the effects of dispersion due to the 
boundary layers, and 2) the standing wave analysis ignores the 
finite impedance of the rigid ends. Of these two differences 
the first is the more important. This results in different 
curvatures of the two theoretical lines. 

The second set of comparisons is made for an open ended 
prime mover used in Che’s thesis research [Ref. 3]. The 
unusual property of this prime mover is that while 1/Q 
decreases with AT for the first and third modes (as is usual 
for a prime mover), 1/Q increases with AT for the second mode. 
This behavior provides a good test for theoretical models. The 
results for 1/0 vs AT are shown in figures 9 through 11 for 
the first, second and third modes, respectively. It is seen 
that the present analysis and the standing wave analysis are 
in close agreement in the first and second modes. The 


agreement diminishes slightly at the third mode. The agreement 


69 


with the measured values is worse than before. However, this 
1s most likely due to contamination of the helium with air and 
water vapor as discuseed in Che’s thesis. Another interesting 
aspect of the analysis is the prediction of an extremum in 1/0 


for the second and third modes. 
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Figure 3. Graph of 1/0 vs AT (K) for the closed end prime 
mover filled with helium at 170 kPa mean pressure. 
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Figure 4. Graph of 1/0 vs AT (K) for the closed end prime 
mover filled with helium at 376 kPa mean pressure. 
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Figure 5. Graph of 1/Q vs AT (K) for the closed end prime 
mover filled with helium at 500 kPa mean pressure. 
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Figure 6. Graph of f, (Hz) vs AT (K) for the closed end prime 
mover filled with helium at 170 kPa mean pressure. 
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Figure 7. Graph of f, (Hz) vs AT (K) for the closed end prime 
mover filled with helium at 376 kPa mean pressure. 
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Figure 8. Graph of f, (Hz) vs AT (K) for the closed end prime 
mover filled with helium at 500 kPa mean pressure. 
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Figure 9. Graph 1/0 vs AT (K) for the first mode of the open 
ended prime mover filled with helium gas at 101 kPa. 
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Figure 10. Graph of 1/Q vs AT (K) for the second mode of the 
open ended prime mover filled with helium gas at 101 kPa. 
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Figure 11. Graph of 1/0 vs AT (K) for the third mode of the 
open ended prime mover filled with helium gas at 101 kPa. 
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Figure 12. Graph of f, (Hz) vs AT (K) for the first mode of 
the open ended prime mover filled with helium gas 


at 101 kPa. 
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Figure 13. Graph of f, (Hz) vs AT (K) for second mode of the 
open ended prime mover filled with helium gas at 101 kPa. 
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the open ended prime mover filled with helium gas 
at 101 kPa. 
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V. SUMMARY 

The results of this thesis can be summarized as follows. 
Arnott’s method was implemented in a Matlab progam. The 
impedance and the pressure distribution were calculated along 
the prime mover. From the frequency response at the left end 
and complex eigenfrequency, the quality factor and the 
resonance frequency was estimated for both cases rigid and 
open end prime movers. The results of the analysis were 
compared with those of a standing wave analysis and with 
experimental results. In the most of the cases the results 
agree well. A future investigation can be the examination of 


the output work of the prime mover. 
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APPENDIX A. THE MATLAB PROGRAM FOR THE OPEN END PRIME MOVER 


% Program for prime mover. 


% PART 1 (CHANGE) 


clear; 

clg; 

coun=0; 6 counter 

pt=0; 6 counter 

qual=zeros (45,1); * quality factor 
inqual=zeros (45,1); $ inverse quality factor 
frO=zeros (45,1); % resonant frequency 
deltaT=zeros (45,1); % temperature difference 
PoOInt=U0; Serna s Cater 

flag=1; % indicator 


slope=sum(’negative’ ); 
load open.dat 6 standing wave analysis data 
$ Main Loop for evrey temperature difference (itrgh) 
fOr TELOGN=293--5-1 50 
coun=1+coun 
numfre=1000; %$ number of frequencies 
6 establish some often used constants 


S—o * number of sections 


J=Sore (—l)e 
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oe 


feu=2/3*(l+eps); 

gamma=5/3* (1l+eps) ; 

spameno- 6. 3143/ (46-3) )* (lteps) 7 

ksolid=.16* (l+eps); 

een =({(1+7) /sqrt (2) ) * (l+eps) ; 

pl=ones (18,numfre) .* (10% (-8)); % presure inside the tube 


fregq=zeros(1l1,numfre) ; % frequncies 


o\e 


w=zeros(1,numfre) ; angular frequencies 


oP 


teren=zeros(1,Ss); end of sections 


oie 


numblay=zeros(1,s); layers of sections 


oe 


lengsec=zeros (1,8); lengths of sections 


tempR=zeros (1,8s); $ rigth temperature of sections 
tempL=zeros (1,8); $ left temperature of sections 
poros=zeros(1,Ss); % porocity of stack 
ratio=zeros(l1,s); % ratio 
SET VALUES 
fremin=864; $ minimum frequency 
fremax=904; $6 maximum frequency 
termin=sum (abs (’free’)),; ¢ end of the last section 
ampres=(1.Q0le+5) *(l+eps) ; $ ambient presure 
dramp=(1.e-16); 6 driver pressure 


section "1" 
terenl=’ opentu’ ; 
teren(1,1)=sum(abs(terenl) ); 


numblay (1,1)=1; 


S5 


lengsec(1,1)=(66.25e-2) * (1l+eps) ; 
tempR (1,1) =293* (1+eps) ; 
tempL (1,1) =293* (1l+eps) ; 
ratio (17d) =(1 2 sige-2) eps, 
poros (1,1) =1* (1+eps) ; 
ares (1, 2) =pi* ratio(l, - je ceps): 
6 section "2" 
teren2='’hexch’; 
teren(1,2)=sum(abs (teren2) ); 
numblay (1,2)=1; 
lengsec (1, 2) =(.82e-2) * (1lt+teps) ; 
tempR (1,2) =293* (1+eps) ; 
tempL (1, 2) =293* (1l+eps) ; 
ratio(1,2)=(5.08e-4) * (l+eps); 
poros (1,2) =.666* (1l+eps) ; 
6 section "3" 
teren3=’ stack’; 
teren(1, 3) =sum(abs (teren3)); 
numblay (1, 3)=11; 
lengsec (1, 3)=(2.59e-2) * (1+eps) ; 
tempR (1, 3) =itrgh* (1+eps) ; 
tempL (1, 3) =293* (l+eps) ; 
ratio(1,3)=(8.636e-4) * (1+eps); 
poros (1, 3) =.89473* (l+eps) ; 


section "4" 


ae 


teren4=' hexch’ ; 
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% 


oe 


teren (1,4) =Ssum(abs (teren4) ); 
numblay (1,4) =1; 

lengsec (1, 4) =(.82e-2) * (l+eps); 
tempR (1, 4) =itrgh* (l+eps) ; 
tempL (1,4) =itrgh* (1+eps) ; 
ratio(1,4)=(5.08e-4) * (1+eps) ; 
poros (1,4) =.666* (l+eps) ; 
section "5" 

teren5=’ opentu’ ; 
teren(1,5)=sum(abs (teren5S5) ); 
mumblay(1,5)=1; 

lengsec (1,5)=(.6852) * (1l+eps) ; 
tempR (1,5) =itrgh* (1l+eps); 
tempL (1,5) =itrgh* (1+eps) ; 
ieee id, 9)— {1.91 7e-2) * (1teps) ; 


poros (1,5)=1* (l+eps) ; 


6 PART 2 
MAIN PROGRAM 


FREQUENCY DISTRIBUTION 


numtot=sum(numblay) +1; 
Pe-pttl, 
tL pt~=1 
adt—=trOteoun=)/ absqual (coun) 


fremin=fr0 (coun-1)-df; 


oy 


fremax=fr0 (coun-1)+df; 
end 
i1=l:numfre; 
if i~=1 
freq=fremin; 
else 
freq(1,i)=fremin+ (fremax-fremin) .*(i./ (numfre-1) ); 
end 
w=2*pi*freq; 
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© GET THE SPECIFIC ACOUSTIC IMPENDANCE AND PRESSURE AT ALL 


\ 


EB Onn. 
% START AT THE RIGHT AN MOVE AT THE LEFT. 
dens (1,numtot) =ampres*4.0e-3 /(tempR(1,s)*8.3143); % density 
visc(1,numtot) =1.887e—5* (tempR(1,s) /27350S)-(. 6567 jae 
viscosity 
kgas (1,numtot) =visc (1,numtot) *cp/npr; % termal condiuicmiyaia, 
sspeed (1,numtot) =972.8*sqrt (tempR(1,s)/273.15); % speed 
%$ isothermal 
typel=’ free’; 
type2=' rigid’; 
if termin==sum (abs (type2) ) 
% IMPEDANCE OF RIGID TERMINATION. 
fac (numtot, 1:numfre)=sqrt ((dens(1,numtot)... 


*sspeed(1,numtot) *2)./(w.*visc(1,numtot))); 
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z(numtot, 1:numfre)=(1+)j) .* (dens(1,numtot) .*sspeed(1,numtot). 
...*fac(numtot, 1:numfre) *sqrt (npr) ) / (sqrt (2) * (gamma-1) ) ; 
% IMPEDANCE OF FREE TERMINATION 
else 
z(numtot, 1:numfre) =free (dens (1,numtot),visc(1,numtot),w, rati 
o(l,s)...,sspeed(1l,numtot),gamma,npr); 
end 
% IMPEDANCE TRANSLATE FOR THE OPEN TUBE OR HEAT EXCHANGER. 
fault=’ stack’; 
truel=’ opentu’ ; 
true2='hexch’ ; 
for k=s:-1:1 
jup=0; 
jJup=numtot-sum(jup+numblay(1,k:s)); 
if teren(1,k)~=sum(abs (fault) ) 
dens (1, jup) =ampres*4.0e-3 /(tempR(1,k).*8.3143); % density 
Vereen, JUup)—1.6e7e—o.* (tempR(1,k)/273.15).%*(.6567); 
+ Viscosiey 
Regen !l, JUup)=visc(1, jup).*cp/npr; % termal conductivity 
sspeed(1, jup) =972.8.*sqrt (tempR(1,k)/273.15) ; %$ speed 
6 isothermal 
S GET LAMBDA AND LAMBDA (T) 
lambda (jup, 1:numfre) =ratio(1,k).*sqrt (dens(1, jup) .*w./visc(1l 
map); 


lambdt (jup, 1:numfre) =sqrt (npr) .*lambda(jup,1:numfre) ; 
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% GET EF(L), EF(L(T)) AND WAVENUMBERS FOR OPEN TUBE 
PARTS 
if teren(1,k)==sum(abs (truel) ) 
flam(jup, 1:numfre) =1- (1+)5) .*sqrt (2) ./lambda ( j up; P2numeeee 
% £(1) 
flamt (jup, 1:numfre) =1-(1+3) .*sqrt (2) ./lambdt (jup, 1:numfre) ; 
pee IE ie) ), 
facl=(1+(gamma-1)/sqrt (npr) )/sqrt (2); 
ka(jup, 1:numfre)=(w./sspeed(1,jup))... 
.* (1+(14)) .*facl./lambda(G@up, tl: numteeee 


else 


oe 


GET F(L), F(L(T)) AND WAVENUMBERS FOR HEAT 

EACHANGERS 1YPE@)  olpiiae 
SQrmi=(ieeg) /soqusie(Z); 
ar(jup,1:numfre) =sqrmi.*lambda (jup, 1:numfre) ; 
argum(jup, 1:numfre) =exp(-ar(jup, 1l:numfre) ); 
ar(jup, 1:numfre) =ar(jup,1:numfre) /2; 

ctanh (jup,1:numfre) =(1l-argum(jup,1:numfre))./(1+argum(jup, 1: 

numfre) ); 

flam(jup, 1:numfre) =l-ctanh(jup,1:numfre) ./ar(jup,1:numfre) ; 
% tis 

ar(jup, 1:numfre) =sqrmi*lambdt (jup,1:numfre) ; 
argum(jup,1:numfre) =exp(-ar(jup,1:numfre) ); 
ar (jup,1:numfre)=ar(jup,1:numfre) /2; 

ctanh (jup, 1:numfre) =(1-argum(jup,1:numfre))./(1+argum(jup,1: 


numfre) ); 
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mame « jUp, l|:numfre)=l—ctanh (jup, 1:numfre) ./ar(jup,1:numfre) ; 
eee Ct }) 
low=w/sspeed(1,  jup); 
ka (jup, 1:numfre) =low. *sqrt ((gamma-(gamma-1)... 
at lamia(gueeeenumtre))./itlam(aup, )snuamire))-; 
end 
% TRANSLATION THEOREM 
comden (jup, 1:numfre) =dens (1, jup)./flam(jup,1:numfre) ; 
Zint (jup, 1:numfre) =comden (jup,1:numfre) .*w./ (ka(jup,1:numfre 
Ipeeeenos (1,k)); 
dsub (1,k)=lengsec(1,k)./numblay (1,k); 
gece, 1 >numire) =ka(jup, 1:numfre) -*dsub (1; k) ; 
sn(jup,1:numfre) =sin(ar(jup,l1:numfre) ); 
cs (jup,1:numfre) =cos (ar (jup,1:numfre) ); 
et (jup, l:numfre) =cs(jup, l:numfre) ./sn(jup,1:numfre) ; 
Paemerguo, l:numfre)=zint (jup,l:numfre) ./z(juptl,1:numfre) ; 
yeu, i snumfre)=zint (gmp, l:-numfre) .* (ct (jup, 1l:numfre)... 
eee s ( jUp,]:numftre)) ./ (facsiqgup, l1:numfre) .*ct (jup, l:numfr 
era): 
pl (jup, 1:numfre) =pl (jup+1,1:numfre).*(cs(jup,1l:numfre)... 
fee acs (jUp, l:numfre) .*sn(jup,1:numfre) ); 
else 
%$ STACK WITH LINEAR DISTRIBUTION OF TEMPERATURE BETWEEN 
TaeeeeNDS OF THE STACK 
toz=(tempR(1,k)-tempL(1,k))/lengsec(1,k); % approximation 


ob gela 
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ns=0; 

for la=ktnumblay (7k) =i: = yk 

ns=nstl; 

tens=tempR(1,k)—-(tempR(1,k)—-tempL(1,k))*ns/numblay (1,k); 
tnsml=tempR(1,k)—-(tempR(1,k)—-tempL(1,k)) * (ns-1) /numblay (1,k) 
tave(l,la)=(tensttnsml) /2; 
dsub(1,k)=lengsec(1,k) /numblay (1,k); 
dens (1,la)=ampres*4.0e-3 /(tave(1,la)*8.3143) ; %$ density 
visc (1,la)=1.887e—-5* (tave(1,la) /273.15)%*(.6567); % vise@euiaey 
kgas (1,la)=visc(1,jup)*cp/npr; % termal conductivity 
sspeed(1,1la)=972.8*sqrt (tave(1,la)/273.15); % speed 

% isothermal 

% GET LAMBDA AND LAMBDA (T) 
lambda (la, 1:numfre) =ratio(l1,k).*sqrt (dens(1,la).*w./visc(1,l 
a)); 
lambdt (la, 1:numfre) =sqrt (npr) *lambda(la,1l:numfre) ; 

ar(la,1:numfre) =sqrmi* lambda (la,1:numfre) ; 

argum(la,1l:numfre)=exp(-ar(la,l:numfre)); 

ar(la,1l:numfre) =ar(la,l:numfre)/2; 
ctanh (la, 1:numfre) =(l-argum(la,1l:numfre))./(l+argum(la,1l:num 
Pee je) y 
flam(la,1l:numfre) =l-ctanh(la,1:numfre)./ar(la,l:numfre) ; 

% £ (1} 

ar(la,1:numfre)=sqrmi*lambdt (la,1:numfre) ; 


argum(la, l1:numfre) =exp(-ar(la, 1:numfre)),; 
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ar(la,1l:numfre) =ar(la,1:numfre) /2; 
ctanh(la,1:numfre) =(l-argum(la,1l:numfre))./(1l+argum(la,1:num 
mee) ) ; 
flamt (la,1l:numfre)=l-ctanh(la,1l:numfre) ./ar(la,1l:numfre) ; 
ee £ (1 (Ee 
low=w/sspeed (1,1a) ; 
ka(la,1:numfre)=low.*sqrt ((gamma-(gamma-1)... 
~.*flamt(la,1l:numfre))./flam(la,1:numfre) ); 
zint (la, 1:numfre) =dens(1,la).*w... 
./ (poros(1,k).*flam(la,1l:numfre) .*ka(la,1l:numfre) ); 
aleocim(la, ];:mumfire) =toz* (flamt (la,l:numfre)... 
7olam(la, l:numfre)-1)./ (2*tave(1,la)*(1l-npr) ); 
kal=ka(la,l:numfre).’; 
Zintl=zint(la,l:numfre).’; 
alprl=alprim(la,1l:numfre).’; 
zeta0=(lengsec (1,k)-(ns-1) *lengsec(1,k)/numblay(1,k)); 
zetafin=(lengsec(1,k)- (ns) *lengsec(1,k)/numblay(1,k)); 
Z0=zZ(latl,l:numfre).’; 
wor—1.e-4; 
[zeta, zout)]=trode45 (’tube’, zeta0, zetafin, z0,tol,kal, zintl,al 
prl); 
n=length (zeta); 
z(la,1l:numfre)=zout(n,:); 
zdumy=z0; 
pl0=pl1 (la+1l,1l:numfre).’; 


tol=1.e-4; 
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{zeta, plout] =trode45(’tubel’, zeta0, zetafin, pl, tol kale 
, ZaOMYy ) 4; 

nl=length (zeta); 

pl(la,1l:numfre) =pHeutinl je), 

end 

end 


end 


+ MINIMIZATION TO GET RESONANCE FREQUENCY AND QUALITY 
BAC TOR: 
tfun=-j*w.*z(1,1:numfre) *dramp./pl1(1,1:numfre); 
fOr a—lias 
pl(a,1l:numfre)=pl (a,1l:numfre) .*tfun; 
end 
amplit=zeros(1,numfre) ; 
amplit (1,1:numfre) =abs (p1(1,1:numfre) ); 
{maxamp, i1)=max(amplit); 
resfre=freq (i) 
amphaf=maxamp/sqrt (2); 
q=0; 
frehaf=0; 
if point== 
for s=2:numfre 
if amplit(s)>=amphaf & amplit (s-1)<=amphaf 
frehaf=(freq(s)+freq(s-1))/2; 


q=0.5*resfre/ (resfre-frehaf) 
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end 
end 
else 
for s=2:numfre 
if amplit(s)<=amphaf & amplit(s-1)>=amphaf 
ama=1 
frehaf=(freq(s)+freq(s-1))/2 
g=0.5*resfre/ (resfre-frehaf) 
end 
end 
end 
if q~=0 
x0=zeros (3,1); 
x0 (1,1) =maxamp; 
x0 (2,1)=resfre; 
x0 (3,1) =q; 
options (1) =0; 
options (2)=1.e-4; 
options (3) =1.e-6; 
options (14) =2000; 
est=fmins(’erro’,x0,options, {],amplit, freq); 
inqual (coun, 1) =est (3,1)%*(-1); 
frO (coun, 1)=est (2,1); 
deltaT (coun, 1) =abs (itrgh-293) ; 
guai(eoun; ))=est (S71); 


end 
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Di Coun-—1 
dif=ingual (coun-l) 1) —-inguatl teounm. ass 
if dif<0 & coun==2 
slope=sum(abs (’positive’)); 
elseif dif<0O & coun==3 
slope=sum(abs (’positive’)); 
end 
1f flag==0 & slope~=sum(abs (’positive’ ) ) 
if dif<=0 | q==0 
peimu—ly 
flag=1; 
for s=2:numfre 
if amplit(s)<=amphaf & amplit (s-—1)>=amphaf 
frehaf=(freq(s)+freq(s-1))/2; 
g=0.5*resfre/ (resfre-frehaf) 
end 
end 
x0 (1,1) =maxamp; 
x0 (2,1)=resfre; 
x0 (3, 1) =a; 
options (1) =0; 
options (2)=1.e-4; 
options (3) =1.e-6; 
options (14) =2000; 
est=fmins'(’ erro’ , x0; optionsmiy amplit, fuecae 


inqual (coun, 1)=est (3,1)%*(-1); 
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fr0(coun,1)=est (2,1); 
deltaT (coun, 1)=abs (itrgh-293) ; 
qual vecun; 1) =est (371) > 
end 
end 


end 


% PART THREE 

%¢ MINIMIZATION USING COMPLEX ANGULAR FREQUENCIES. 

eptions (1) =0; 

options (2)=1.e-4; 

options (3)=1.e-4; 

options (14) =700; 

we 2 er te (Coun, 1) *(1—7/ (2 Guat eoun, §))) ; 
y2=wlt+wl1*10.e-6; 
y l=wl-wl1*10.e-6; 

ella=fmin (’erro3’,yl,y2, options, itrgh) ; 
fr0 (coun, 1) =real(ella)/(2*pi); 

qual (coun, 1)=-real (ella) /(2*imag(ella)); 
inqual (coun, 1) =qual (coun, 1)%*(-1); 


end 


ae 


PLOTS AND STORE THE DATA (CHANGE). 
v—ieenmgqual(iscounys)) 7 (delta; coun, 1) )’; 
eae cCoun,1))’ ); 

Eldg—pepen(’ ingfis.uxt’,’wt’); 


Mmerintet (fid,”’ Ingual=%—-9.6f£ T=%-4.1f RF=%-5.2f\n’',y); 
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leoad-conpeZ dat 


plot (deltaT(1:coun) , inqual (1:coun) ,compf2 (271) compi& Aree 


figure 


ao 


, Open (pf eweopen(: 7,4). 4") 
title(’Inverse Quality Factor V.S Temperature 
Difference.Open End’); 
epertrel 
xlabel(’Negative Delta "T" ’); 
ylabel(’ 1/0"); 
gtext (‘Ambient pressure 101 Kpa’) 
gtext (‘Experimental Result "o"’) 
gtext (‘Computer Result "*"’) 
gtext (‘Third Mode’ ) 
print -deps inguf3 
figure 
plot (delvar(Tveoun) £ru @: coun) 
xlabel(’Negative Delta "T" ’); 
ylabel(’Resonance Frequence’ ); 
title(’Resonance Frequency V.S Temperature 
Difference.Open End’); 
gtext (‘Ambient pressure 101 Kpa’) 
gtext (’Third Mode’ ) 
@nasrel 


print -deps frf3 
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f 


% FUNCTION TRODE45.M 


PuNnct ion 


feomernyout )=trode45 (FunFcn,t0,tfinal,y0,tol,kal,zintl,alprl) 


pepema—(1/4 3/8 12/13 1 1/2)’; 


beta=[ [ i 0 0 0 
[ 5 9 0 0 
(igeZ2 =7Z00 9296 0 


P8341 —326 sz 29440r= 845 


0 


[G08 0 41040263525 97955 = 5043 


0)/4 
Oy 32 
O] 7/2197 
0] /4104 


Oy 20S20 |; 


gamma=[ (902880 0 3953664 3855735 -1371249 277020] /7618050 


[Pr=Z210 9070 22528 21970 
hae 

Peet lalization 
ete 0; 
hmax= (t-tfinal) /5; 
himan= (t-tfinal) 7/2000; 

hn (e—-t final) / 100; 
y=y0(:); 
f=y*zeros (1,6); 
Seut—t; 
yout=y./’; 
tau=tol*max(norm(y,’inf’),1); 
$ Main loop 
while (t>tfinal) & (n>=hmin) 


Pewee ternial! ~h—t—-tiinal-. end 
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= 0 484275 60)/ 752400 


% Compute the slopes 
temp=feval (Funken, kal, zinti,alpriy 
£(:,1)=temp(:); 

EOL j=1:5 
temp=feval (FunFcn, kal, zintl, alprl,t-alpha(j)*h, y-h*f*beta(:, 
De 
£(:, j+1)=temp(:); 
end 
$ Estimate the error and the acceptaple error 
delta=norm(nh*f*gamma(:,2),’inf’); 
Eau=tol*max (norm (yy int’); 1)? 
% Update the solution only if the error is acceptaple 
if delta<=tau 
tt hi, 
y=y—-h*f*gamma (:,1) ; 
joie = eels 216) 5 
VOuE=(yout sy. |) 
end 

% Update the step size 

pow=1/5; 
if delta~=0.0 

h=min (hmax,0.8*h* (tau/delta) *pow) ; 
end 

end 
VE (et inet) 


disp( ’singularity likely.’) 
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ee = 


le 

end 

% FUNCTION "TUBE.M" GIVES THE DERIVATIVE OF IMPEDANCE 
function dzdz=tube (kal, zintl,alprl, zeta, z) 


eee Kale *zintl.* (1—(z./zintl) .*2)+2.*alprilg*z); 


% FUNCTION "TUBE1.M" GIVES THE DERIVATIVE OF PRESSURE 
function dpdz=tubel (kal, zintl, zdumy, zeta,p1) 
epaz—).*pl.*kal.*zintl./zdumy; 


Bonkoeb or OR RMINS . LEAST SQUARS TECHNIQUE. 
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See UNCT LON 


function rmserr=erro(x0,amplit, freq) 


mamp=x0 (1,1); 

EV—KO(Z, 1); 

Z=x0 (3,1); 
focl=(mamp*f0) ./ (z*freq); 


pmmemrnboc)].°~2)./((f£0./freq-freq/f£0) .*2+1/2°2) )-amplit.%2; 


rmserr=sum(err.%*2); 


2oRROS.M. Okan MIN 2s LEAST SOUARS TECHNIQUE 
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*¢ FUNCTION 
function fplus=erro3 (wl, itrgh) 


[z,dens, sspeed, visc]=argo3f (wl,itrgh); 
gamma=5/3* (l+eps) ; 
npr=2/3* (l+eps); 


See dAdo sort (2) )*(ik+teps); 


De 


° 


— 


oO 


[z 


fofw=sqri*sqrt (dens*3*sspeed“4*npr... 
./(wl.*visc))./(gamma-1) ; 

sl=real (fofw) ; 

S2=real (z); 

s3=imag (fofw) ; 

s4=imag(z) ; 


fplus=(Ss1">2—-s24270me a(S 3°2-S4-2). 2, 


FUNCTION “ARGOSE JMS USE DS Biaer Reo. MoM 


funcevon 
, dens, sspeed, visc])=argo3 (wl, itrgh) 
w=wl; 


establish some often used constants 


s=5; $ number of sections 
J=Scibe (1) 
npr=2/3* (lt+eps) ; 
gamma=5/3* (l+eps) ; 
Gp=(279* 8 99M 3/"(4 Fee) ) *OErepsn 
ksolid=.16* (1l+eps) ; 
sqri=((14+4) /sqre (Zp (1 tepse 
teren=zeros(s,1); 
numblay=zeros(1,Ss); 


lengsec=zeros (1,S); 


BZ 


tempR=zeros (1,S); 
tempL=zeros(1,Ss); 
poros=zeros(1,8s); 
ratio=zeros (1,8); 

SET VALUES 
termin=sum (abs (’ free’)); 
ampres=(1.0le+5) * (l+teps) ; 


dramp=(1.e-8) ; 


Section “i" 

terenl=’ opentu’ ; 
teren(1,1)=sum(abs (terenl)); 
numblay (1,1)=1; 

lengsec (1,1)=(65.34e-2) * (l+eps) ; 
tempR (1,1) =293* (1l+eps) ; 

tempL (1,1) =293* (1+eps) ; 
ratio(1,1)=(1.917e-2) * (1l+eps) ; 
poros (1,1) =1* (1+eps) ; 
ares(1,:)=pi*ratio(1l,:)*(l+eps); 
Section "2" 

teren2=' hexch’ ; 
teren(1,2)=sum(abs (teren2)); 
MmumoLay (i, 2)=1; 

lengsec (1,2) =(.82e-2) * (lt+teps) ; 


tempR (1,2) =293* (l+eps) ; 


BS 


tempL (1, 2) =293* (1+eps) ; 
ratio(1,2)=(5.08e-4) * (1l+eps) ; 
poros (172) —. 650674 ttepair, 

% SeCCEIOnN 3. 
teren3=' stack’; 
teren (1,3) =sum(abs (teren3)); 
numblay (1, 3)=11; 
lengsec (1,3)=(2.59e-2) * (1+eps) ; 
tempR (1, 3) =itrgh* (1+eps) ; 
tempL (1,3) =293* (1+eps) ; 
ratio(1, 3)=(8.636e-4) * (1l+eps); 
poros (1,3) =.89473* (1l+eps) ; 

6 section "4" 
teren4=’ hexch’ ; 
teren(1,4)=sum(abs (teren4)); 
numblay (1, 4)=1; 
lengsec (1,4) —(.82e=2) 7 (iteps)® 
tempR (1,4) =itrgh* (1l+eps) ; 
tempL (1,4) =itrgh* (l+eps) ; 
ratio(1,4)=(5.08e-4) * (1l+eps); 
poros (1, 4) =.666* (1l+eps) ; 

6 section "5" 
teren5=' opentu’; 
teren(1,5)=sum(abs (teren5) ); 
numblay (1,5) =1; 


lengsec(1,5)=(.6852) *(l+eps); 
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tempR(1,5)=itrgh* (1+eps) ; 
tempL (1,5) =itrgh* (1+eps) ; 
mateo 105) ='( 1.91 7e—-2) (itepsi; 
poros (1,5) =1* (l+eps) ; 
1f s>=7 
6 section "6" 
teren6=’hexch’; 
teren(1,6)=sum(abs (teren6)); 
numblay (1, 6)=1; 
lengsec(1,6)=(1.02e-2) * (1+eps) ; 
tempR (1, 6) =293* (1+eps) ; 
tempL (1, 6) =293* (1+eps) ; 
pamo(1,6)=(1.02e-3) * (1l+eps) >; 
poros (1, 6) =.667* (1+eps) ; 
Sesecmion “7" 
teren7=' opentu’; 
teren(l1, 7)=sum(abs (teren?7) ); 
numblay (1, 7)=1; 
lengsec (1,7)=(5.483e-2) * (1l+eps) ; 
tempR (1,7) =293* (1+eps) ; 
tempL (1, 7) =293* (1l+eps) ; 
ratio(1,7)=1.917e-2* (lt+eps) ; 
poros(1,7)=1* (lt+eps) ; 


end 


* START AT THE RIGHT AND MOVE AT THE LEFT. 


SD 


numtot=sum(numblay) +1; 
dens (1, numtot) =ampres*4.0e-3 / (tempR (1, s) *8.3143); % density 
visc(1,numtot) =1.887e-5* (tempR(1,s) /273.15) * (.6567) ; 
$ vViseoscaay 
kgas (1, numtot) =visc (1,numtot) *cp/npr; % termal conductivity 
sspeed(1,numtot) =972.8*sqrt (tempR(1,s)/273.15); %* speed 
%$ isothermal 
typel=’ free’; 
type2=' rigid’ ; 
1f termin==sum (abs (type2) ) 
% IMPEDANCE OF RIGID TERMINATION. 
fac3(1,numtot) =sqrt ( (dens{1, numtom ere 
*sspeed(1,numtot) .*2)./(w. *visc (1,numtot) )); 
z(1,numtot)=(1+35) .* (dens (1,numtot).*sspeed(1,numtot)... 
.*fac3(1,numtot) *sqrt (npr) ) / (sqrt (2) * (gamma-1) ) ; 
6 IMPEDANCE OF FREE TERMINATION 
else 
z(1,numtot) =free (dens (1, numtot),visc(1,numtot),w, ratio(1,s). 
,sspeed(1l,numtot), gamma, npr) ; 


end 


%¢ IMPEDANCE TRANSLATE FOR THE OPEN TUBE OR HEAT EXCHANGER. 
bante=' stack, 
truel=’ opentu’ ; 
true2=' hexch’ ; 


for k=s:-1:]1 
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jup=0; 
jJup=numtot-sum(jupt+numblay (1,k:s)); 
if teren(1,k) ~=sum(abs (fault) ) 
dens (1, jup) =ampres*4.0e-3 /(tempR(1,k).*8.3143); % density 
wirce, 10p)=1.887e-5.* (tempR(1,k) /273.15) .*(.6567) ; 
Se VISCOSIEY 

kgas(1, jup)=visc(1,jup).*cp/npr; % termal conductivity 
sspeed(1, jup) =972.8.*sqrt (tempR(1,k)/273.15); % speed 


% isothermal 


% GET LAMBDA AND LAMBDA (T) 
Pmeada(l, Jup)=ratio(1,k).*sqrt (dens (1, jJup) .*w./visc(1, jup)); 
lambdt (1, jup) =sqrt (npr) .* lambda (1, jup) ; 
% GET F(L), F(L(T)) AND WAVENUMBERS FOR OPEN TUBE 
PARTS 
1f teren(1,k)==sum(abs (truel) ) 
meaner, yup) =1—(1+ 3) .*sqrt (2) ./lambda (1, jup) ; % ene} 
meanerul, JUuD)=1—-(1+ )) .*sqrt (2) ./lambdt(1,jup); 3 £(1(t)) 
fac31=(1+ (gamma-1)/sqrt (npr) )/sqrt (2) ; 
ka(1,jup)=(w./sspeed(1,jup))... 
Pate + ye bacol 27 Lambada (1, jup) ) ; 
else 
% GET F(L), F(L(T)) AND WAVENUMBERS FOR HEAT 
EXCHANGERS TYPE "SLIT". 
sgqrmi=(1- 4) /sqrt (2); 


ar (1, jup) =sqrmi.*lambda (1, jup) ; 


| 


argum (1, gup) —exp (-ar (iy jum oe. 
ar (1, jup)=ar (1, jup) /2; 
ctanh (1, jup) =(l-argum(1, jup)) ./ (1l+argum@ gues, 
flam(1, jup) =l-ctanh (1, jup)e/ar @sy iup); + ufe( 1) 
ar(1, jup) =sqrmi*lambdadt (1, jup) ; 
argum(1, jup) =exp(-ar(1,jup)); 
ar(1,Jjup)=ar(1, jup) /2; 
etanh (1, jup) =(1-argum(1, jup)) ./ (1+argum (1a 
flamt (1, jup) =1-ctanh (1, jup)./ar(1,jup); See 
low=w/sspeed(1, jup); 
ka(1l, jup)=low.*sqrt ((gamma-(gamma-1)... 
.*flamt (1, jup))./flam (ifm 
end 
6 TRANSLATION THEOREM 
comden (1, jup) =dens (1, jup) ./flam(1, jup) ; 
zint (1, jup) =comden (1, jup) .*w./ (ka(1, jup) .*pores(i7 ae 
dsub (1,k) =lengsec(1,k)./numblay(1,k); 
ar(1,jup) =ka(1, jup) .*dsub(1,k); 
sn(1, jup)=sin(ar(1,jup)); 
cs(1, jup)=cos(ar(1,jup)); 
et (1, jup) =cs(1, jupye/snG july 
fac3 (1, jup) =zZint@, jup) 4/72 awe 
Z(1, jup) =zint (lagu) aaiict (17 jue) eae 
-j.*fac3(1, jup)) ./(fac3 (1, jup) -*ct (1, ¢upieeee 


else 
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% STACK WITH LINEAR DISTRIBUTION BETWEEN THE ENDS OF THE 
DEACK. 
toz=(tempR(1,k)—-tempL(1,k))/lengsec (1,k); $ approx 
ay At 
ns=0; 
fom lLa—ktnumblay (1,k)-1:-1:k 
ns=nstl; 
tens=tempR(1,k)-(tempR(1,k)-tempL(1,k))*ns/numblay (1,k) ; 
tnsml=tempR(1,k)-(tempR(1,k)-tempL(1,k)) * (ns-1) /numblay (1,k) 
tave(1,la)=(tens+tnsml) /2; 
dsub (1,k)=lengsec (1,k) /numblay (1,k); 
dens (1,la)=ampres*4.0e-3 /(tave(1,la)*8.3143); $ density 
Vise (1,la)=1.887e-5* (tave(1,la)/273.15)%*(.6567); 
viscosity 
kgas (1, la)=visc(1,jup)*cp/npr; % termal conductivity 
sspeed(1,la)=972.8*sqrt (tave(1,la)/273.15); $ speed 
6 isothermal 
% GET LAMBDA AND LAMBDA (T) 
lamia(l, la)=ratio(1,k).*sqrt(dens(1,la) .*w./visc(1,la)); 
lambdt (1, 1a) =sqrt (npr) *lambda(1,1a); 
ar(1,la)=sqrmi*lambda(1,1a); 
argum(1,la)=exp(-ar(1,la)); 
ar(1,la)=ar(1,la)/2; 
Cctanh(1,la)=(l-argum(1,la)) ./(lt+targum(1,1la)); 


flam(1,la)=l-ctanh(1,la)./ar(1,la); soe ele) 


ao 


ar(1,la)=sqrmi*lambdt (1,1la); 
argum(1,la)=exp(-ar(1,la)); 
ar(1,la)=ar(1,1la)/2; 
ctanh(1,1la)=(l-argum(1,la))./(1+argum(1,la)); 
flamt (1, la)=l-ctanh (1, tape ar (1, 1a); % £ (1 (ee 
low=w/sspeed(1,1a); 

ka(1,la)=low.*sqrt ((gamma-(gamma-1)... 

.*flamt(1,la))./flam(1,la)); 
Zint (1,la)=dens(1,la).*w... 
./ (poros (1,k) .*flam(1,la) .*ka( tages 
alprim(1,la)=toz* (flamt (1,1a) 
./flam(1,la)—-1) ./(2*tave (1, la)* (1-npaee 
Kal-kai(l, lajeee, 
Zintl=Zaiiee (ly ae 
alprl=alprim(l1,la).’; 
zeta0=(lengsec(1,k)—-(ns-1) *lengsec(1,k) /numblay(1,k)); 
zetafin=(lengsec(1,k)-(ns) *lengsec (1,k) /numblay (1,k)); 
zO=z(1,lt+la).’; 
tol=1.e-4; 
[zeta, zout]=trode45 (’tube’, zeta0, zetafin, z0,tol,kal,zintl,al 
prl); 
n=length (zout) ; 
Z41, la) =zoul(n ae 
end 
end 


end 
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ZZ tl i) ys 
dens=dens (1,1); 
sspeed=sspeed (1,1); 


visc=visc(l1,1); 
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